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Molecular simulations were carried out to evaluate the performance of intermolecular potential models for the 
adsorbates: nitrogen, ethylene and carbon dioxide with leading quadrupoles. Models having only LJ interactions, 
are compared with the corresponding multi-site models with partial charges for their performance in describing 
the experimental data for vapour-liquid equilibrium (VLE), adsorption isotherms and isosteric heat of adsorption 
on planar adsorbent substrates. It is found that the models with partial charges outperform models using only LJ 
atoms for the description of adsorption, highlighting the important role of the quadrupole moment in a non- 
uniform adsorbed phase. These models were then used to construct the wetting maps of nitrogen, ethylene 


and carbon dioxide as a function of the temperature and the strength of the planar substrate to show the in- 
fluence of the quadrupole moment on the behaviour of non-wetting/partial wetting/complete wetting at the 
coexistence pressure of the bulk phase. 





1. Introduction 


Adsorption of gases on graphitized thermal carbon black (GTCB) has 
been extensively studied because this solid possesses an energetically 
homogeneous surface and is ideal to characterize the intrinsic interac- 
tion between adsorbate and the surface. They also have relevance to the 
adsorption behaviour of carbon adsorbents used in practice. Adsorbates 
can form either molecular layers or clusters on the surface, depending on 
the interplay of the various potential energy interactions in the system. 
To understand the growth mode of the adsorbed film or the behaviour of 
clusters on the surface, extensive studies were carried out to investigate 
the phenomenon of wetting/non-wetting on GTCB, for example Ar 
[1-2], Kr [3-4], Xe [5], quadrupolar adsorbates, Nz [6], and C2H4 [7], 
COz [8], and polar adsorbates such as H20 [9-10], NH3 [11], CoCloF4 
[12], CHgCl, [13]. A number of theories on wetting/non-wetting have 
been proposed in the literature [14-19], and because of the different 
terminologies used in different scientific disciplines to describe the same 
phenomena we provide brief definitions of the various terms in Table 1 
to give clarity to the subsequent discussion. It is known that wetting/ 
non-wetting is governed by the interplay between the adsorbate inter- 
molecular interactions, the adsorbate-adsorbent interactions and the 
temperature [17,20]. However, the role of molecular dipole moments 
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and quadrupole moments on wetting/non-wetting behaviour in 
adsorption has not been sufficiently quantified. Terlain and Larher [21] 
reported isotherms for CO2, N20 and C2Nz2 adsorption on graphite, and 
found the quadrupole moment strongly favours an ordering orientation 
within the adsorbate, and as such a question is raised about its role in 
wetting/non-wetting phenomena. 

Fig. 1 compares the strength of the quadrupoles for the adsorbates 
chosen in this investigation with argon as the reference. For these ad- 
sorbates, we chose a number of potential models: Single LJ-site, Two-LJ- 
site, and Multi-Site with fixed partial charges to illustrate not only the 
role of the quadrupole moment but also the molecular shape on the 
adsorption/wetting behaviour. Computer simulations of vapour liquid 
equilibria and adsorption of Nz, C2H4 and CO% on graphite have been 
carried out to evaluate these models against the experimental data and 
simulations for the adsorption of these gases on planar substrates having 
different affinities. The results have been used to construct wetting maps 
and compare these with maps for non-polar argon. 
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Table 1 
Definitions of the various terms used. 
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Term Definition 





Bulk coexistence pressure The bulk coexistence pressure (Po) refers to the 
saturation pressure for vapour-liquid equilibrium at 
temperatures above the bulk triple point 
temperature, or the sublimation pressure for vapour- 
solid equilibrium at temperatures below the bulk 
triple point temperature. 

The adsorbate sparsely covers the surface at the bulk 
coexistence pressure (Po). Typically, the adsorbed 
density at Po is less than the statistical monolayer 
density. 

A system is called incomplete wetting when the 
adsorbate forms a film of finite thickness at Po. This 
is called Type 2 [24] or Stranski-Krastanov growth 
[25] in the colloid science literature. 

A system is called complete wetting when the 
adsorbed density tends to infinity when the pressure 
approaches Po. This is called Type 1 [24] or Frank- 
van der Merwe growth [25]. 

The temperature at which there is a first order 
transition in the adsorbed density at Po. If this 
transition leads to complete wetting (i.e. the 
adsorbed density tends to infinity) the temperature 
is denoted as the wetting temperature Tw. If a new 
layer n, is formed on top of an adsorbed film the 
temperature of layer formation is called the layering 
temperature of the n-th layer, Ti, n- 

The temperature, denoted as Tr, which demarcates 
the wetting region from the incomplete wetting 
region. Above Tx ,the interface between the 
adsorbed film and the gas phase becomes 
sufficiently undulated such that complete wetting 
occurs, otherwise incomplete wetting occurs at 
temperatures less than Tp. 


Non-wetting 


Incomplete wetting 


Complete wetting 


Wetting temperature/ 


layering temperature 


Roughening temperature 





Ar N, C,H, co, 
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Fig. 1. The absolute value of the quadrupole moment (10° C m°) for Ar(0), 
CO2(14.98) [22], No(4.70) [23] and C2H,4(6.67) [23]. 


2. Theory 
2.1. Kinetic Monte Carlo 


The kinetic Monte Carlo (kMC) method was used in this study 
because of its distinct advantage of determining the chemical potential 
accurately [26], and has been successfully used to describe vapour- 
liquid equilibria [27] and adsorption [26,28]. Simulations were car- 
ried out with the substrate positioned at the bottom of a rectangular 
simulation box and a hard wall at the top of the box. Periodic boundary 
conditions were applied in the x- and y- directions, in the plane of the 
substrate. 10° configurations were sampled in each of the equilibration 
and production stages. 


2.2. Grand Monte Carlo simulation 


Grand ensemble Monte Carlo simulation (GCMC) with the Metrop- 
olis scheme was used in simulations to derive the wetting map with one 
billion configurations each in the equilibration stage and in the pro- 
duction stage. A new configuration in the Markov chain, was selected 
from displacement, insertion or deletion with equal probability. In the 
equilibration stage, the maximum displacement length was initially set 
as half of the smallest dimension of the box and adjusted to achieve an 
acceptance ratio of 20%, and then kept constant in the production stage. 

The surface excess, the isosteric heat, the orientation distribution and 
the local density distribution obtained from the simulation are defined 











Table 2 
Molecular Parameters of adsorbates used in this study. 
Model Atom o(nm) eg /ke(K)  q(e) Geometry 
Ar M-—Ar Ar 0.3405 119.8 0 - 
[29] 
N2 2LJ + 3Q N 0.331 36.0 —0.482 (N —N) = 
[31] COM* 0.0 0.0 0.964 0.110nm 
1CLJ [30] - 0.3615 101.5 0 0 
C2H4 2LJ + 4Q CH2 0.3575 99.8 0.32 (C—C) = 
[33] 0.170nm 
= - - —0.32  I(q-q) = 
0.130nm 
2LJ [32] CH2 0.3675 85 0 (C—C) = 
0.133nm 
CO2 3LJ + 3Q Cc 0.28 27.0 0.7 C—O) = 
[31] o 0.305 79.0 —0.35 0.116nm 
1CLJ [35] - 0.36481 246.15 0 -= 





* COM = centre of mass; **: the partial charges are on the line perpendicular to 
the C-C bond. 


in Appendix A. 
2.3. Fluid-Fluid (Upp) and Solid-Fluid (Usp) interaction energies 


The molecular parameters for the potential models used to describe 
the Ar, N2, C2H4, and COz interactions, are listed in Table 2. The model 
for argon is taken from Michels [29], and denoted as M—Ar. Nitrogen 
was modelled either as a single LJ molecule (1CLJ N2) [30] or as two-LJ 
sites with three fixed partial charges (2LJ + 3Q Ng) [31], Ethylene was 
modelled as a two-site LJ model (2LJ C2H4) [32] or by the two-LJ sites 
and four partial charges (2LJ + 4Q C2H4) [33]. The single LJ-site model 
of C2H4 was not chosen here because it fails to describe the experimental 
vapour-liquid equilibria (VLE) [34]. CO2 is modelled with either a single 
LJ site model (1CLJ CO2) [35] or with three LJ-sites plus three partial 
charges (3LJ + 3Q CO2) [31]. 

We use the ratio of the minimum of the solid—fluid (usr) potential 
energy profile to the minimum of the fluid—fluid (urp) potential energy 
profile, D*, to define the relative strength of an adsorbent with respect to 
a given adsorbate: 








(1) 


The intermolecular energy was calculated as the sum of all pairwise 
interaction energies with the cross molecular parameters computed 
from the Lorentz-Berthelot mixing rules. The minimum fluid-fluid po- 
tential energy of a pair of adsorbate molecules is calculated as a function 
of the separation distance between their centres of mass. For a given 
separation, both molecules were given 10° random orientations, and the 
minimum energy from the set of pairwise potential energies was 
recorded for this separation. These energies were then plotted against 
separation distance, and the minimum of this potential profile, Urr min, 
used in eq. (1). 

The solid-fluid potential energy between an LJ site and a planar 
graphite substrate, modelled as a continuum solid, was calculated from 
the 10—4-3 potential equation [36], as a function of the separation dis- 
tance in the z- direction between the centre of mass of a given adsorbate 
and the solid surface. 

The adsorption strength of a substrate is denoted by the parameter 


D= |usr-ynin / | ure min 


Dss(K) = p,02,€ss, which only depends on the molecular properties of the 
substrate. For graphite, D,, = 124 K using the standard values €,, = 28 K 
and o = 0.34 nm and p = 38.2 nm 2. To model planar substrates of 
different surface affinity, we varied the well-depth of the interaction 
energy, &ss/kp. 


3. Results and discussion 


Simulations with the various intermolecular potentials listed in 
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Fig. 2. Simulated isotherms with 1CLJ and 2LJ + 3Q for Nə adsorption at 60 K (a) and 77 K (b) on graphite. Experimental data [37] is shown as the black solid line. 


Values of Po are given in Table 8. 





Fig. 3. Snapshots of molecular configurations for nitrogen with (a) 1CLJ model and (b) 2LJ + 3Q model at 60 K and 700 Pa. The red circle encloses a 


“pinwheel” cluster. 


Table 2 generally give good agreement with experiment for the vapor- 
-liquid equilibria (Fig. 17, Fig. 18, Fig. 19), with the models that account 
for shape and quadrupole moment performing slightly better. We next 
consider their performance in the description of isotherm and isosteric 
heat for adsorption on graphite. 


3.1. Adsorption on graphite 


N2 

The simulated isotherms and the experimental data [37] for N2 
adsorption on graphite at 60 K and 77 K are shown in Fig. 2 to illustrate 
the difference in the adsorption behaviour below and above the bulk 
triple point temperature of 63 K. 

The 2LJ + 3Q N2 model with partial charges describes the experi- 
mental data at 60 K and 77 K better than the 1CLJ model which is 
particularly poor at 60 K, where the molecular shape and quadrupole 
moment play a significant part in determining the packing of the 
structured first layer and consequently have a propagation effect in 
building the higher layers. The 1CLJ model over-predicts the adsorbed 
density at higher loadings because the 1CLJ molecules pack more effi- 
ciently than linear molecules as a CP 2D layer. The configuration 
snapshots of the monolayer at 60 K in Fig. 3 illustrate the differences in 
packing for the two models. 

The regular hexagonal pattern of the 1CLJ molecules, provides 


adsorption sites at the centres of three adjacent molecules, on which 
molecules adsorb to form the second layer. This regular makes the sec- 
ond layer adsorption sites stronger, and as a result the 1CLJ model over- 
predicts the adsorbed density, compared to the experimental data. On 
the other hand, the 2LJ + 3Q model provides a more irregular mono- 
layer structure, with pinwheel formations [38-39] in which a central 
molecule is oriented normal to the surface and its six nearest molecules 
lie flat on the surface (for example the cluster located inside the circle in 
Fig. 3b), as noted earlier by Kuchta and and co-workers [39]. 

Fig. 4 shows the isosteric heat versus loading at 60 K, obtained from 
the simulations with the 2LJ + 3Q model. The simulation results capture 
the pattern of the experimental data in the sub-monolayer region, but 
there are deviations between the simulation results and the experi- 
mental data, which we attribute to the difference between the simula- 
tion results at 60 K and experimental results for nitrogen adsorption on 
highly graphitized CarbopackF [37]. 


3.1.1. C2H4 

Ethylene has a higher quadrupole moment than nitrogen and we 
therefore expect to see a stronger effect due to the presence of the 
quadrupole. The simulated isotherm for ethylene adsorption on graphite 
at 105 K (the bulk triple point) and the experimental data of Menaucourt 
et al. [40] are shown in Fig. 5. 

To improve agreement between experiment and simulation in the 
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Fig. 4. Simulated isosteric heat with the 2LJ + 3Q Nz model (triangle symbols) 
and experimental data [37] (the solid line) is for nitrogen adsorption at 60 K 
on graphite. 


Henry law region, the cross-molecular parameter between the adsorbate 
and the adsorbent, computed by the Lorentz-Berthelot rule, was modi- 
fied to include a binary interaction parameter ky [41], eş = 


(1 —ksf),/ (esep), with k = 0.15. 


The simulated isotherms over the multilayer coverage range for the 
2LJ + 4Q model are in better agreement with the experimental data, 
than those for the LJ model, which supports the argument that the 
quadrupole moment needs to be included in the description of a non- 
uniform adsorbed film and affects the adsorbate packing in higher 
layers. 

The distribution of molecules for the 2LJ + 4Q potential model, at P/ 
Po = 0.1, at which the monolayer is completed, is shown in Fig. 5b as a 
3D-plot of density versus distance from the graphite surface and the 
angle between the molecular axis and the normal vector. The monolayer 


(a) s0 C,H;-105K 





2LJ+4Q 
ksf=0.15 


Surface excess (umol/m?) 








0.0 0.2 0.4 0.6 0.8 1.0 


PIP, 


Chemical Engineering Journal 419 (2021) 129502 


is about 0.4 nm from the surface and the 90 degree peak in the distri- 
bution confirms that most molecules lie flat on the surface with the 
remaining molecules adopting a vertical orientation, similar to that 
observed for nitrogen. 

Fig. 6a shows the simulated isosteric heats versus loading with the 
2LJ + 4Q model together with the component contributions from the 
adsorbate-adsorbent (SF) interactions and the adsorbate—adsorbate (FF) 
interactions. 

Within the sub-monolayer coverage region (adsorbed density less 
than 8 pmol/m2), the isosteric heat is greater than the heat of conden- 
sation, and is constant across the first-order transition in the first layer 
which has been discussed in an earlier paper [42]. The importance of the 
electrostatic terms can be seen in Fig. 6b where the adsorbate—adsorbate 
contribution is decomposed into LJ and electrostatic contributions. 
Before the first-order condensation of the first layer occurs, the contri- 
bution from the electrostatic interactions to the heat released by the 
adsorbate—adsorbate interactions is less than 2 kJ/mol, and it increases 
very sharply to 6 kJ/mol at the completion of the monolayer, as the 
molecules rearrange and re-orient in the monolayer in order to maxi- 
mize the quadrupole—-quadrupole interactions, with some molecules 
adopting the vertical orientation seen in the orientation distribution in 
Fig. 5b. 


3.1.2. CO2 

Simulations using the 1CLJ and 3LJ + 3Q models adsorbed on 
graphite at 193 K and 273 K, were tested against the data of Beebe et al. 
[43] and of Guillot and Stoeckli [44]. The 1CLJ model over-predicts the 
saturation pressure [45] at 193 K, but the simulated adsorption isotherm 
agrees with the experimental data if plotted against absolute pressure as 
shown in the inset in Fig. 7a. The simulated isotherms and the experi- 
mental data at 193 K show that the adsorbed density at the coexistence 
pressure Po is less than twice the statistical monolayer density (incom- 
plete wetting), while those at 273 K, plotted in Fig. 7b, show complete 
wetting. 

Fig. 8 shows the simulated isosteric heats obtained with the 1CLJ and 
3LJ + 3Q models together with the experimental data of Spencer et al 
[43] at 193 K. The 3LJ + 3Q CO2 model fits the experimental data better 
because the adsorbate—adsorbate interactions (FF) are augmented by the 
electrostatic contributions which optimise the molecular packing at the 
completion of the monolayer (around 8 pmol/m?). 

From the analysis of the isotherm and the isosteric heat for N2, C2H4, 
and CO, adsorption on graphite we can draw the following general 
conclusions: 


(b) 2LJ+4Q C,H, 














Local Density(-) 





90 


30 o 
0 angie ) 


Fig. 5. (a) Simulated isotherms with 2LJ and 2LJ + 4Q models for ethylene adsorption at 105 K on graphite. Experimental data is shown as the black solid line. The 
value of Po is given in Table 9. (b) Orientation density distribution of the centre of mass at 105 K and 0.1 reduced pressure. 
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Fig. 6. (a) Simulated isosteric heat and contributions from the adsorbate-adsorbent (SF) and adsorbate—adsorbate interactions (FF). (b) The contributions from (LJ) 
and electrostatic interactions (CL) to the adsorbate-adsorbate contributions (FF) for C2H4 adsorption on graphite at 105 K. 
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Fig. 7. Simulated isotherms with 1CLJ and 3LJ + 3Q models for CO» adsorption at (a) 193 K and (b) 273 K on graphite with kf = 0.05 and 0.0635 for the 1CLJ model 
and the 3LJ + 3Q model [41 ]respectively. The inset figure is the isotherm plotted versus absolute pressure. Experimental data are shown as the black solid line. The 


value of Po is given in Table 10. 


1. The experimental isotherm and the isosteric heat are better fitted 
when electrostatic interactions (quadrupolar models) are included in 
the potential functions. 

2. The packing of molecules in the first layer is strongly affected by the 
quadrupole moment, which has a propagation effect in building 
higher layers. 

3. Most molecules in the first layer on the surface lie flat on the surface, 
with a few molecules adopting other orientations; for example in the 
pinwheel structures of nitrogen. 


3.2. The effects of adsorbate quadrupoles on wetting maps 


The quadrupolar potential models for nitrogen, ethylene and carbon 
dioxide have been used to simulate adsorption on planar substrates of 
different adsorption strength to construct wetting maps. The choice of 
saturation pressure, Po, is critical in the construction of a wetting map, 
and the details of our Pg calculations for the adsorbates in this paper are 
given in Table 7 to Table 10. 


3.2.1. Generic mechanism of wetting behaviour 

Using GCMC simulation, we obtained wetting/non-wetting maps, 
and explored the microscopic origins behind the layering temperature 
and the wetting temperature for argon, nitrogen, ethylene and carbon 
dioxide. The common features, shared by all these adsorbates, are shown 
in the schematic wetting and heat maps in Fig. 9a and 9b, respectively. 
From these maps we distinguish two adsorbent substrates to exemplify a 
strong adsorbent I, and a weak adsorbent HI, located by the yellow circles 
in Fig. 9. 

The wetting map in Fig. 9a shows regions of non-wetting, incomplete 
wetting and complete wetting over the space bounded by temperature 
and adsorbent strength. These regions are separated by three lines: (1) 
the wetting temperature (Ty) line (blue dashed line) separating the 
wetting and non-wetting regions, (2) the roughening temperature (Tr) 
line (the red dotted line) separating the wetting and incomplete wetting 
regions and (3) the layering temperature line for the first layer (Tz1) 
(black solid line); separating the incomplete wetting and non-wetting 
regions. Within the incomplete wetting region, there is a family of 
layering temperatures for the second and higher layers, with Tr n+ > 
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Fig. 8. Simulated isosteric heats and the contributions from the adsorbate-adsorbent (SF) and adsorbate—adsorbate (FF) interactions; (a) 1CLJ model and (b) 3LJ + 


3Q model for CO2 adsorbed on graphite at 193 K. 
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Fig. 9. (a) Schematic map showing the separation of non-wetting/incomplete wetting/complete wetting zones in a temperature and surface affinity Dss space. (b) 
The schematic map for the heat of sublimation (condensation) of fluid molecule and the isosteric heats at zero loading as a function of temperature. 


Tın Where n is the layer n. From this wetting map, the layering tem- 
perature of the first layer for adsorbent I, T;1, is 0 K, with stronger 
substrates having T;,; = OKand weaker adsorbents having T,,1 > OK. For 
substrate II the layering temperature becomes the wetting temperature, 
which is the same as the roughening temperature. 

The schematic map in Fig. 9b shows the heats of sublimation (Asub) 
and condensation (Acond) at temperatures above and below the triple 
point temperature of the bulk phase, and the isosteric heats at zero 


loading (q®) for adsorbates on substrates of various strength as a 
function of temperature. As the loading is increased the adsorbate 
configuration contributes to the subsequent adsorption of higher 
adsorbate layers and therefore plays a part in the wetting/non-wetting 
behaviour at Po. 

To elucidate the distinction between substrates I and II, we establish 
the connection between the wetting map and the heat map for these 
substrates. We define the difference between the isosteric heat at zero 


loading qo and the heat of sublimation or condensation at 
temperatureT: 


Asl T) = q® — As A cona (T) = 4) — Acona 


When either Asub Or Acond is positive, incomplete wetting occurs and 
transition to wetting occurs when the temperature is greater than the 
roughening temperature. 

We denote the energy contributed from the adsorbate—adsorbate 
interactions during the course of adsorption by Qer. 

Substrate I has a layering temperature T;,; of 0 K in the first layer. 
Between T;,; and Tr there is incomplete wetting, and complete wetting 
occurs at temperatures above the roughening temperature. This sub- 
strate is at the lower limit of strong adsorbents, and any weaker sub- 
strates have a layering temperature greater than 0 K in the first layer, 
which increases with decreasing affinity. 

Fig. 9b shows that at temperatures below the triple point tempera- 
ture A,,,»(T)is less than 0 because the isosteric heat at zero loading is less 
than the heat of sublimation. However, the contribution Qpr, from the 
adsorbate—adsorbate interactions is >A,,,(T), and therefore the layering 
temperature for the first layer is 0 K. As the temperature is increased, the 
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Fig. 10. (a) The wetting map in T-D,, space for Ar showing (b) Plots of the heat of sublimation (solid circles) for Ar obtained from simulations of vapour-solid 
equilibria by Chen et al. [45], the heat of condensation (empty circles) obtained with kMC simulations of vapour-liquid equilibria. [47]. The solid lines show the 
isosteric heats at zero loading obtained by Monte Carlo integration of the Boltzmann factor [48] for substrates I and II as marked in the wetting map (a). 


difference between Asu(T) and Qrr becomes smaller and higher layers 
are formed sequentially in the incomplete wetting region. When the 
temperature has reached the triple point temperature (the roughening 
temperature), complete wetting occurs because the probability of 
adsorbate separations close to the potential minimum increases as the 
amplitude of thermal fluctuations increases. 

For substrates having a higher affinity than substrate I, the layering 
temperature of the first layer Tı, is 0 K because the magnitude of 


Asub(0) = qo —Asyp is smaller than Qpr, as shown for substrate III in 
Fig. 9a. Argon, nitrogen and ethylene adsorption on graphite are ex- 
amples of Substrate II. 

On the other hand, for substrates having weaker affinity than sub- 
strate I, the layering temperature for the first layer T,,; is greater than 0 
K because A,,(0)is negative and is not exceeded by Qpr. Therefore the 
first layer does not form unless the temperature is higher as is the case 
for substrate IV in Fig. 9a. Carbon dioxide adsorbed on graphite is an 
example of substrate IV. 

Substrate II is an example in which the layering temperature of the 
first layer is the same as the wetting temperature at the roughening 
temperature. When the first layer is formed, higher layers are also 
formed because of the high thermal fluctuations; and Trı = Tw = 
Tr = Ti. The heat maps add further insight to understanding the con- 
ditions controlling wetting. Plots of Q against T show that when the 
temperature has reached the triple point temperature, adsorption on the 
surface competes with the cohesiveness of the bulk liquid, and that the 


difference Agona(Tir) = qo —A-ona then becomes equal to or greater than 
Qrr the energy contributed by the adsorbate—adsorbate interactions, 
resulting in a complete wetting. 

For substrates having affinity less than that of substrate II, the wet- 
ting temperature is greater than the roughening temperature because 
the adsorbent strength is too weak to induce wetting unless there is 
sufficient thermal fluctuation. An example is substrate V. 

We next consider the actual wetting maps for the adsorbates exam- 
ined in this paper. The values of D,,(K) for substrates I and II are sum- 
marized in Table 11. 


3.2.1.1. Argon. The non-wetting/wetting behaviour of adsorbed argon 
on planar graphite substrates has been previously studied in detail by 
Prasetyo et al. [46] using the M—Ar potential model, and their results 
are summarized here in the wetting map in Fig. 10a. The substrates I and 
II discussed above are also labelled in the wetting and heat maps in 


Fig. 10a and Fig. 10b. 

To further investigate the role of Qrr in the formation of the first 
layer, we simulated the adsorption of argon adsorption on 10-4-3 
graphite at 40 K. At this temperature the isotherm exhibits a first-order 
transition. Fig. 11a and 11b show the GCMC and kMC isotherms and the 
isosteric heat curve for argon adsorption on graphite. 

Snapshots from the GCMC isotherm during the build-up of the first 
layer, are shown in Fig. 11c for points A to E marked on the isotherm. 
The layering process starts with the formation of clusters of molecules 
(Point A), followed by the coalescence of clusters to form circular 2D- 
islands (Point B). As more molecules are added the 2D-island becomes 
a 2D-dense phase with (approximately) linear interfaces separating it 
from the rarefied phase (point C). Further adsorption occurs at the linear 
boundaries between the two phases. At this point, Qpr is due to the 
interaction between the added molecule and molecules at the boundary 
of the condensed phase, and has a constant value of, Qrr = 2.2kJ/mol, 
shown as a red curve in Fig. 11b. 

As adsorption progresses and the number of molecules in the 
condensed phase become larger, the rarefied phase becomes a mono- 
layer surrounding a 2D-vapour bubble (point D). Finally, the monolayer 
is completed as a hexagonal packing, as confirmed by the 2D-radial 
distribution at E (Fig. 11d). The mechanism described here for molec- 
ular layering across a vdW loop of the canonical isotherm is applicable 
for the other adsorption systems considered in this paper. 

We can now establish a connection between the wetting map and the 
heat map to explain the mechanisms behind the substrates I and II 


1. Substrate I: At 0 K, complete wetting occurs. 

2. Substrate IJ: At temperatures between 0 K and the triple point tem- 
perature, the difference A,,»(T) is negative and the system is non- 
wetting. At the triple point temperature the adsorbate—adsorbate 
interactions, reach a value of Qrr = 2.2kJ/mol which is equal to 
-Asub(T) , and complete wetting occurs. 


3.2.1.2. N2. The wetting map and heat map from the simulations of N2 
adsorption on planar substrates of various affinities are summarized in 
Fig. 12a and b. The differences Asu(T)and Acona(T), at O K and the triple 
point temperature for the Substrates I and HI are listed in Table 3. The 
energy provided by the adsorbate—adsorbate interactions prior to the 
formation of the first layer, Qrr, as 2.6 kJ/mol. Details of this are pre- 
sented in Fig. 21. Table 3 shows that this energy is > -Asat 0 K for 
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Fig. 11. (a) Isotherms for argon adsorption at 40 K on graphite, simulated in GCMC and kMC-NVT ensembles. (b) Isosteric heat and separate contributions from the 
adsorbate—adsorbate interactions (FF) and the adsorbate-adsorbent interactions (SF) calculated by canonical ensemble simulation [49]. (c) Snapshots of molecules in 
the first layer at the points marked on the canonical isotherm. (d) 2D radial distribution in the first layer (the inset is the snapshot of molecules for point E). 


Substrate I and therefore the layering temperature of the first layer for 
this substrate is 0 K. For Substrate II, the energy from the adsorbate- 
-adsorbate interactions balances A,,,q at the triple point temperature, 
and hence Tı = Tw = Tr = Tp for Substrate II. 


3.2.1.3. C2H4. The wetting map and the heat maps for ethylene 
adsorbed on 10—4-3 graphite are shown in Fig. 13a and b, respectively. 
Substrates I to IT are marked in these figures and the details of the 
isosteric heat versus loading for ethylene adsorption at 50 K are shown 
in Fig. 21. The energy contributed from the adsorbate—adsorbate in- 
teractions in building the first layer, Qrp, is 6.3 kJ/mol and was found to 
be >-Asu (0) for Substrate I and -A;onaat the triple point temperature for 
Substrate II given in Table 4. 


3.2.1.4. CO2. Fig. 14a and Fig. 14b show the non-wetting/wetting and 
heat maps for CO% adsorption on substrates of different strength with the 
locations of Substrates I to I marked. The values of A,,5(0),Acona(Tir) in 
Table 5 show that the adsorbate—adsorbate energy is > Asub or Acona for 
substrate I and therefore Tı = OK. For Substrate II, T,1 = Tw = Tr = 
Ti since Qrris > Acona(Tir). The decomposed heat curves are shown in 
Fig. 21. 

Xu et al. [8], have shown that the bulk triple point temperature 
separates the incomplete wetting region from the complete wetting 


region. 


3.2.2. Comparison of wetting maps for adsorbate molecules 

We now consider the effects of electrostatic interactions on the 
wetting behaviour of the adsorbates on substrates of different strengths. 
Values of D* (Gr), defined in Eq. (1), are given in Table 6 and show that 
D* (Gr) decreases with increase in the quadrupole strength , suggesting 
that the strength of the adsorbate-adsorbent interaction is lower as 
quadrupole strength increases. compared to the intermolecular in- 
teractions. For argon and nitrogen, graphite is a stronger adsorbent than 
Substrate II because the isosteric heat at zero loading is higher than the 
heat of sublimation, therefore, the layering temperature for the first 
layer is 0 K. Ethylene has a stronger quadrupole moment than nitrogen, 
its affinity for graphite is higher than Substrate II; therefore the layering 
temperature of the first layer is O K even though the isosteric heat at zero 
loading is less than the heat of sublimation. Here the difference is made 
up from the energy provided by the adsorbate—adsorbate interactions. 
Incomplete wetting occurs for temperatures between 0 K and the 
roughening temperature of 105 K, which is in agreement with the 
experimental data of Sutton et al. [51] and Menaucourt et al. [40]. CO2 
which has the strongest quadrupole moment of the four adsorbates, 
graphite falls in the class of moderately strong adsorbents, i.e. between 
Substrate I and Substrate II and the layering temperature for the first 
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Fig. 12. (a) The wetting map for 2LJ + 3Q N2 showing different wetting regions as a function of the temperature and the surface affinity Dss. (b) Plots of the heat of 
sublimation (solid circles) for 2LJ + 3Q N2 from Brown and Ziegler [50], the heat of condensation (empty circles) are from kMC simulations of vapour-liquid 
equilibria [47]. The solid lines show the isosteric heats at zero loading for substrates of different strength obtained by Monte Carlo integration of the Boltzmann 


factor [48]. 


Table 3 

The difference between the heat of sublimation or condensation and isosteric 
heat at zero loading for Nz on chosen substrates (I and II) and its maximum 
adsorbate—adsorbate interaction energy. 








Substrate Asup(0) (kJ/mol) Acond (Tir) (kJ/mol) Qer (kJ/mol) 
I —2.2 —0.5 2.6 
II —3.8 —2.6 


layer for graphite was found to be 90 K. 

To compare the wetting behaviour of adsorbates that have different 
intermolecular strengths, as reflected in the triple point and critical 
temperatures, we define a non-dimensional temperature: 


0 = (T—T,)/(Te — Tr) 
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The triple point and critical temperatures for these adsorbates are 
given in Table 6. The scaled roughening temperature, Or = (Tr — T,r)/ 
(T. — Ty) are also shown in the same table. 

The parameter D* (eq. (1)) is the ratio of the strength of the 
adsorbate-adsorbent interactions to that of the adsorbate—adsorbate 
interactions for a given adsorbent. To assess the adsorbent strength of a 
substrate, we take graphite as a reference and define a reduced affinity 


Table 4 

The difference between the heat of sublimation or condensation and isosteric 
heat at zero loading for C2H4 on chosen substrates (I and II) and the maximum 
adsorbate—adsorbate energy. 








Substrate Asub(0)(kJ/mol) Acona(Tir) (kJ/mol) Qrr(kJ/mol) 
I —4.3 —0.24 6.3 
II —10.0 —6.20 
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Fig. 13. (a) The wetting map for ethylene as a function of temperature and adsorbent strength Dss. (b) Plots of the heats of sublimation (solid circles) for 2LJ + 4Q 
C2H, from Brown and Ziegler [50], and the heats of condensation (empty circle) obtained from kMC simulations of VLE [47]. The solid lines show the isosteric heats 
at zero loading for substrates of different strength obtained by Monte Carlo integration of the Boltzmann factor [48]. 
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Fig. 14. (a) The wetting map for carbon dioxide as a function of the temperature and the adsorbent strength Dss. (b) Plots of the heat of sublimation (solid circles) for 
3LJ + 3Q CO: obtained from simulations of vapour-solid equilibria by Chen et al. [45], the heat of condensation (empty circles) obtained from kMC simulations of 
vapour-liquid equilibria [47]. The solid lines show the isosteric heats at zero loading for substrates of different strength obtained by Monte Carlo integration of the 


Boltzmann factor [48]. 


Table 5 

The difference between the heat of sublimation or condensation and isosteric 
heat at zero loading for CO2 on chosen substrates (I and II) and its maximum 
adsorbate—adsorbate interaction energy. 











Substrate Asub (0) (kJ/mol) Acona( Tir) (kJ/mol) Qer (kJ/mol) 
I —4.3 6.12 8.1 
I —18.5 —7.8 

Table 6 


Physical properties of the linear adsorbates (the units of the quadrupole moment 
are 10°° Cm’). 




















Fluid Quadrupole Tir/(K) Tc/(K) D*(Gr) Tr/(K) Or 
moment 
Ar 0 82.5 154.5 9.0 80 —0.04 
[45] [45] 
N2 —4.7 63 126.5 8.5 66 0.05 
[31] 
C2H4 6.7 104 286.3 4.6 105 0.01 
[32] 
CO2 —12.2 212 306.8 3.3 215 0.03 
[45] [45] 
Table 7 
Summary of Po of M—Ar at different temperature used in our 
simulations. 
Temperature (K) Po(Pa)for M—Ar 
40 2.04 x 10° 
Table 8 


Summary of Po of nitrogen for different temperature used in simulations. 


= D 


D* (Gr) 

The wetting maps reconstructed in a reduced temperature 8 reduced 
affinity space are shown in Fig. 15. 

For a given reduced temperature 0; as marked in Fig. 15, we see that 
non-wetting is expected for all adsorbates on very weak substrates (less 
than A). As the affinity of the substrate increases complete wetting oc- 
curs for argon and nitrogen for substrates with affinities between A and 
B but ethylene and carbon dioxide are non-wetting, due to the effects of 


Table 9 
Summary of Po for ethylene for temperature used in simulations. 





Temperature (K) Po(Pa)for Exp Po(Pa)for 2LJ Po(Pa)for 2LJ + 4Q 





105 1.46 x 10? 4.86 x 102 1.45 x 10? 





Table 10 
Summary of Po of carbon dioxide (3LJ + 3Q model) for the different tempera- 
tures used in this paper. 





Temperature (K) Po(Pa)for Exp Po(Pa)for 1CLJ Po(Pa)for 3LJ + 3Q 








193 9.67 x 104 3.72 x 10° 1.06 x 105 
273 3.48 x 10° 3.53 x 10° 3.50 x 10° 
Table 11 


Strength of the substrate (D,,(K)for Ar, N2, C2H4 and COz on substrates I and II. 








Substrate adsorbate molecule 

Ar No C2H4 CO2 
I 55.8 42.4 100.3 227.95 
I 21.3 16.6 39.7 32.54 








Temperature (K) Po(Pa)for Exp Po(Pa)for 1CLJ Po(Pa)for 2LJ + 3Q 
60 5.94 x 10° 6.40 x 10° 6.20 x 10° 
77 9.31 x 10+ 9.13 x 104 9.50 x 104 
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Fig. 15. Wetting maps of argon, nitrogen, ethylene and carbon dioxide in terms of the reduced temperature and Q. 
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Fig. 16. (a) Local density distribution of the centre of mass of M—Ar at 83 K, 2LJ + 3Q N2 at 77 K, 2LJ + 4Q C2H4, at 116 K and 3LJ + 3Q CO at 215 K at Po versus 


distance away from the graphite surface. 3D plots of orientation density distributions of (b) 2LJ + 3Q N2 at 77 K, (c) 2LJ + 4Q C2H4, at 116 K and (d) 3LJ + 3Q CO2 at 
215 K at Po on a graphite surface. 
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Fig. 17. (a) The vapour-liquid coexistence curve for N2. (b) Logarithm of the saturated 
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vapour pressure versus the inverse of temperature for N2. Experimental data 


are shown as the solid line and the simulation results from the 2LJ + 3Q Nz model and the 1CLJ Nz model are shown as triangles and circles, respectively. 
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Fig. 18. (a) The vapour-liquid coexistence curve for C2H4. (b) Logarithm of the saturated vapour pressure versus the inverse of temperature for C2H4. Experimental 


data are shown as solid lines and the simulation results for the 2LJ + 4Q C2H4 model 
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and the 2LJ C2H4 model are shown as the triangles and circles, respectively. 
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Fig. 19. (a) The vapour-liquid coexistence curve for CO2. (b) Logarithm of the saturated vapour pressure versus the inverse of temperature for CO2. Experimental 
data are shown as solid lines and the simulation results from the 3LJ + 3Q CO2 model and the 1CLJ CO2 model are shown as triangles and circles, respectively. 
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Fig. 20. Logarithm of the sublimation pressure versus the inverse of temperature for Nə. The values calculated from empirical expressions for sublimation pressure 
[50] are shown as a dashed line with triangles and the simulation results for the 2LJ + 3Q Nz model with Bin-Monte Carlo [54] are shown as a solid line with circles. 
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Fig. 21. Isosteric heat of adsorption of (a)argon at 40 K, (b) nitrogen at 20 K,(c) ethylene at 50 K and carbon dioxide at 90 K on graphite surface and its contributions 
from the Fluid- Fluid (FF) and Solid-Fluid (SF) interactions to the first layer. The plateau of the fluid—fluid contribution is the adsorbate—adsorbate interaction energy 


provided to build the first layer,Qrr. 
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the quadrupoles as discussed earlier. When the affinity is greater than C, 
wetting occurs for all adsorbates. 

The local density distributions for argon, nitrogen, ethylene and 
carbon dioxide adsorption on graphite at their respective bulk coexis- 
tence pressures and roughening temperatures are plotted in Fig. 16a as a 
function of the distance between the solid and the centre of the fluid 
molecule. For argon, the peaks are sharp, with clear separation between 
molecular layers. On the other hand, the first peaks for nitrogen, 
ethylene and carbon dioxide are broader with a shoulder, indicating of 
some molecules have vertical orientations. The 3D-orientation density 
distributions plots for nitrogen, ethylene and carbon dioxide in Fig. 16b, 
c, d, respectively confirm this structure. Most molecules adopt orienta- 
tions parallel to the surface to maximize their LJ interactions with 
graphite. This orientation with respect to the surface is common for all 
linear adsorbates on graphite at low temperatures [8,37,52] and in 
carbon pores [34], as confirmed with the analysis of heat of adsorption 
by Hoory and Prausnitz [53]. 


4. Conclusions 


The simulations reported in this paper have address three issues, 
relating to the wetting on a planar substrates. First, simulations with 
potential models accounting for the presence of a quadrupole moment 


Appendix A 


The surface excess density in the system, T is defined as 


(N ) — VaccP 


S 


r= 
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give a good account of both vapour-liquid equilibria and adsorption. 
Secondly, models which include electrostatic interactions have the same 
characteristics of non-wetting, incomplete wetting and complete wet- 
ting, and are governed by the difference between the isosteric heat at 
zero loading for adsorption on planar substrates and the heats of subli- 
mation and condensation of the bulk fluid. Finally, it was found that 
most quadrupolar molecules adopt orientations parallel to the surface to 
maximize their interactions with graphite. 
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where (N) is the ensemble average of the number particles in the system, Vaccis the accessible volume, where solid-fluid potential is non-positive, 


p,is the density of the bulk gas,S is the surface area of the solid. 
The isosteric heat, qst, is calculated from the following equation : 


(No) f(U,N) 
FNN) © FN, N) = (Vace/Vo)f(No,No) 





qst = 


where N, is the number of molecules occupying the accessible volume of the system at the same density as the bulk phase, N is the number of 
molecules in the system, Vacc is the accessible volume within which the solid-fluid potential energy is non-positive, Vp is the volume that used to 
simulate the bulk gas phase, U is the configurational energy of the system,f (X, Y) is the fluctuation variable, defined as f(X, Y) = (X Y} —(X)(Y), and (-) 


denotes an ensemble average. 


The orientation of a fluid molecule is defined as the angle, 8, between the linear molecular axis and the normal vector through the graphite surface. 


The orientational density distribution was calculated as: 


(Nz) 
a chy A 
POO = FFA sind) AG 


where (N,,) is the average number of molecules in the region bounded by [z,z + Az] and angles are bounded by [0,0 + A4]. 
The local density in the distance of the centre of mas of adsorbate molecules from the surface was calculated as: 


= (AN, 2442) 
p(z) = STAS 


where (N,z+az) is the average number of molecules whose centre of mass in the region bounded by [z,z + Az]. 


6 Appendix B. Physical properties 


6.1. Vapour-liquid equilibrium 


Fig. 17 to Fig. 19 show the performance of the various potential models for the description of the vapour-liquid equilibria (VLE) as plots of the 
temperature versus the coexistence densities and the vapour pressure as a function of temperature for N2, C2H4 and CO2. These simulation results were 
obtained with kMC simulations. It is seen that the multi-site models with fixed partial charges yield a better description of the phase diagram and the 
vapour pressure than the models with only LJ potentials, indicating that molecular shape and the electrostatic contributions to the potential are 


significant. 


14 


H. Xuet al. Chemical Engineering Journal 419 (2021) 129502 


6.2. Bulk coexistence pressure 


6.2.1. Ar 
For argon, the bulk coexistence pressure, Pg was determined from kMC vapour-liquid simulations [47]. At temperatures above and below the triple 
point temperature, we used the simulation results for vapour-solid equilibrium by Chen et al. [45]. 


6.2.2. No 

For the 2LJ + 3Q potential model for N potential, the bulk coexistence pressure, Pg, was determined from kMC simulations of the vapour-liquid 
equilibria [47] at temperatures above the bulk triple point temperature, Tr. Values of Pg at temperatures below Tr, were taken from experimental data 
[50] since these are in excellent agreement with the simulation results of vapour-solid equilibria [54] as confirmed by the plots given in Fig. 20. 


6.2.3. CoH4 

With the 2LJ + 4Q potential model for C2H4, we obtained the coexistence pressure of the bulk phase, Po, from kMC simulations of vapour-liquid 
equilibria (VLE) [47] at temperatures greater than the triple point temperature. , Below the triple point temperature we used the results from Smukala 
et al. [55]. The experimental data are taken from the NIST [56]. 


6.2.4. CO2 

The coexistence pressure of the bulk phase, Pg, was determined from kMC simulations of vapour-liquid equilibrium [47] at temperatures above the 
triple point temperature. Below the triple point temperature we used the simulation results for vapour-solid equilibrium from Chen et al. [45]. The 
experimental data were taken from the NIST [56]. 
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